Contexto
Esta serie consta del valor FOB en dólares de las importaciones que llegan a los puertos de Colombia vía marítima y su destino final es la ciudad de Bogotá, en el periodo transcurrido entre enero del 2012 hasta diciembre de 2021, la serie es de tipo mensual.
Definición: El valor FOB en dólares de la mercancía, es valor de la mercancía en el momento que se carga a bordo del medio de transporte, en este caso el marítimo.
La serie consta de 120 observaciones, lo que corresponde a los 10 transcurridos desde el 2012 hasta el 2021
# 1.000'000.000
vafodo <- ts(importaciones[,3], start = c(2012, 01), frequency =12)/1000000000
plot(vafodo, ylab = "Miles de millones de dólares", main = "Valor FOB", lw =2)
Visualmente vemos que la serie presenta una tendencia, la cual parece ser creciente con el tiempo. El rango de valores que toma la variable se va haciendo cada vez mayor (heteroscedasticidad).
1. Parte descriptiva
1.1 Estabilización de la varianza
Transformación de Box-Cox
serie <- vafodo
a <- MASS::boxcox(lm(serie ~ 1), seq(-1, 1, length = 50))
BC.m <- a$x[which.max(a$y)]
BC.f <- forecast::BoxCox.lambda(serie, method = "loglik",
lower = -1,
upper = 1)
# Transformación logarítmica
lserie <- log(vafodo)
a <- MASS::boxcox(lm(lserie ~ 1), seq(-2, 2, length = 50))
BC.ml <- a$x[which.max(a$y)]
BC.fl <- forecast::BoxCox.lambda(lserie, method = "loglik",
lower = -2,
upper = 2)
c(BC.f, BC.fl, BC.m, BC.ml)## [1] -0.2500000 -0.2000000 -0.1919192 -0.1414141Los valores de \(\lambda\) obtenidos por el método de Box-Cox tanto en el paquete MASS, como en el paquete forecast, son diferentes de 1 e inferiorese a 0, además el intervalo de confianza para \(\lambda\) no captura el 1, por lo que se usará la transformación logarítmica. Una vez aplicada, notamos que el intervalo de confianza caputra al 1, por otra parte los valores de \(\lambda\) siguen siendo inferiores a 0.
1.2 Estimación de la tendencia
En el primer gráfico podemos ver que el IC para \(\lambda\) no captura al 1, además el valor de \(\lambda\) para estabilizar la varianza es -0.25, por tanto usaremos \(\lambda =0\) para estabilizar la varianza. En el segundo gráfico ya podemos ver esta transformación logarítmica aplicada, ahora en el IC está incluido el 1 y el \(\lambda\) que estabiliza la varianza es -0.2.
1.2.1 Estimación de la tendencia (lineal)
fit_lserie <- lm(lserie ~time(lserie), na.action = NULL)
summary(fit_lserie)##
## Call:
## lm(formula = lserie ~ time(lserie), na.action = NULL)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.6748 -0.1490 0.0355 0.2354 0.5725
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.251e+02 1.920e+01 -16.93 <2e-16 ***
## time(lserie) 1.630e-01 9.521e-03 17.12 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3011 on 118 degrees of freedom
## Multiple R-squared: 0.713, Adjusted R-squared: 0.7106
## F-statistic: 293.2 on 1 and 118 DF, p-value: < 2.2e-16# Grafico
plot(lserie, lw = 2, main = "Valor FOB en escala log", ylab = "Log de miles de millones de dólares")
abline(fit_lserie, col = "blue", lw = 2)
# Eliminando la tendencia
lserie.sin.tend <- lserie - predict(fit_lserie)
plot(lserie.sin.tend, main = "Valor FOB en escala log sin tendencia", lw =2)
acf(lserie, lag.max = length(lserie))
acf(lserie.sin.tend, lag.max = length(lserie.sin.tend))
La estimación de la tendencia por medio de una regresión lineal simple nos da los parámetros \(\hat{\beta}_0\) y \(\hat{\beta}_1\) significativos, sin embargo al observar la gráfica tenemos que una regresión lineal simple no es una forma óptima de eliminar la tendencia.
1.2.2 Promedio móvil
descomposicion_lserie <- decompose(lserie)
plot(descomposicion_lserie)
1.2.3 Tendencia desde el STL
Como se vio en el punto anterior los datos no presentan una tendencia del todo lineal, además no se veían indicios de estacionalidad, por lo que es necesario aplicar el filtro de promedio móvil. Al aplicarlo podemos ver que la tendencia definitivamente no es lineal, la componente residual no muestra un comportamiento estacional.
En el primer gráfico se hizo un ajuste STL sin ajustar los parámetros, para el segundo se ajustó un polinomio de grado 2, el cual se ajusta mejor que la anterior.
indice_lserie <- as.Date(as.yearmon(tk_index(lserie)))
indice_lserie1 <- yearmonth(as.yearmon(tk_index(lserie)))
# Forma alternativa de extraer el indice
df_lserie <- data.frame(Fecha = indice_lserie,
lserie = as.matrix(lserie))
tibble_lserie <- tibble(df_lserie)
tsibble_lserie <- as_tsibble(df_lserie)
# Primera aproximación al ajuste STL
tsibble_lserie %>%
timetk::plot_time_series(Fecha, lserie,
.interactive = TRUE,
.plotly_slider = TRUE)#<<<<<<< HEAD
# Ajuste STL
# escala log
tibble_lserie %>%
mutate(lserie_ajust = smooth_vec(lserie,
span = 0.3,
degree = 2))## # A tibble: 120 × 3
## Fecha lserie lserie_ajust
## <date> <dbl> <dbl>
## 1 2012-01-01 2.76 2.90
## 2 2012-02-01 2.88 2.92
## 3 2012-03-01 2.98 2.94
## 4 2012-04-01 2.77 2.95
## 5 2012-05-01 3.16 2.97
## 6 2012-06-01 3.16 2.99
## 7 2012-07-01 3.13 3.00
## 8 2012-08-01 3.18 3.02
## 9 2012-09-01 3.06 3.04
## 10 2012-10-01 2.99 3.05
## # ℹ 110 more rows#=======
# Ajuste STL
tibble_lserie %>%
mutate(lserie_ajust = smooth_vec(lserie,
span = 0.2,
#>>>>>>> main
degree =2)
)## # A tibble: 120 × 3
## Fecha lserie lserie_ajust
## <date> <dbl> <dbl>
## 1 2012-01-01 2.76 2.86
## 2 2012-02-01 2.88 2.90
## 3 2012-03-01 2.98 2.93
## 4 2012-04-01 2.77 2.96
## 5 2012-05-01 3.16 2.99
## 6 2012-06-01 3.16 3.00
## 7 2012-07-01 3.13 3.02
## 8 2012-08-01 3.18 3.04
## 9 2012-09-01 3.06 3.04
## 10 2012-10-01 2.99 3.04
## # ℹ 110 more rows# Ajuste STL moviendo los parámetros
# escala log
tsibble_lserie %>% mutate(
lserie_ajus = smooth_vec(lserie, span = 0.3, degree = 2)) %>%
ggplot(aes(Fecha, lserie)) +
geom_line(size =1.05)+
geom_line(aes(y = lserie_ajus), color = "blue", size =1.05) +
theme_bw()
tsibble_lserie %>% mutate(
lserie_ajus = smooth_vec(lserie, span = 0.2, degree = 2)) %>%
ggplot(aes(Fecha, lserie)) +
geom_line(size =1.05)+
geom_line(aes(y = lserie_ajus), color = "blue", size =1.05) +
theme_bw()
tsibble_lserie %>% mutate(
lserie_ajus = smooth_vec(lserie, span = 0.2, degree = 2),
dlserie_ajus = lserie - lserie_ajus) %>%
ggplot(aes(Fecha, dlserie_ajus)) +
geom_line(size =1.05)+
theme_bw()
1.2.4 STL Tendencia y estacionalidad
tsibble_lserie <- as_tsibble(lserie)
tsibble_lserie %>%
model(
STL(value ~ trend() +
season(window = "periodic"),
robust = TRUE)) %>%
components() %>%
autoplot() +
theme_minimal()
1.3 Diferencia Ordinaria
# escala log
tsibble_lserie|>mutate(
diff_lserie = tsibble::difference(value, lag = 1,
differences = 1))|>
autoplot(.vars = diff_lserie, size = 1.05) +
labs(subtitle = "Cambios en escala log del valor FOB") +
theme_bw()
tsibble_lserie <- tsibble_lserie|>mutate(
diff_lserie = tsibble::difference(value, lag = 1,
difference = 1))
# Diferenciando basado en el objeto tibble
tibble_lserie %>%
mutate(diff_lserie = lserie - lag(lserie)) %>%
plot_time_series(Fecha, diff_lserie)tibble_lserie <- tibble_lserie %>%
mutate(diff_lserie = lserie - lag(lserie))
dlserie <- diff(lserie)1.4 Relaciones no lineales dispersión
par(mar = c(3,2,3,2))
astsa::lag1.plot(dlserie, 12, corr = T)
1.5 ACF
acf(dlserie, 48, main = "Serie diferenciada y con logaritmo del valor FOB")
pacf(dlserie, 48)
acf(dlserie, lag.max = 50, main = "Serie diferenciada y con logaritmo del valor FOB")
El acf no parece sugerir una componente estacional. ## 1.6 Índice
AMI
par(mar = c(3,2,3,2))
astsa::lag1.plot(lserie, 12, corr = F)
nonlinearTseries::mutualInformation(lserie, lag.max = 100,
n.partitions = 50,
units = "Bits",
do.plot = TRUE)
## $time.lag
## [1] 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
## [19] 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35
## [37] 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53
## [55] 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71
## [73] 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89
## [91] 90 91 92 93 94 95 96 97 98 99 100
##
## $mutual.information
## [1] 4.862322 3.052912 2.910639 2.853770 2.920942 2.891607 2.907651 2.827147
## [9] 2.925675 2.864164 2.945057 2.907872 2.930455 2.937367 2.914099 2.921245
## [17] 2.897458 2.858685 2.728650 2.655587 2.641206 2.651833 2.831493 2.653982
## [25] 2.687045 2.755120 2.669964 2.698785 2.654698 2.684328 2.620457 2.686898
## [33] 2.690152 2.653327 2.639095 2.577666 2.568432 2.634660 2.650707 2.454465
## [41] 2.437976 2.415092 2.475064 2.416767 2.337317 2.276132 2.396371 2.311431
## [49] 2.352713 2.306940 2.231579 2.096361 2.185558 2.277735 2.168338 2.190601
## [57] 2.096567 2.011887 1.924846 1.933693 1.876562 1.929006 1.939626 1.783764
## [65] 1.850553 2.089983 1.879139 1.926220 1.970048 1.892796 1.973020 1.991408
## [73] 2.021388 1.989870 1.978963 2.251709 2.071896 2.081328 2.122268 2.166044
## [81] 2.207717 2.342629 2.320075 2.513545 2.468303 2.236186 2.380430 2.596091
## [89] 2.451598 2.517136 2.321928 2.476648 2.340975 2.253077 2.179329 2.402292
## [97] 2.312907 2.239678 2.344462 2.271873 2.395462
##
## $units
## [1] "Bits"
##
## $n.partitions
## [1] 50
##
## attr(,"class")
## [1] "mutualInf"# sobre la serie diferenciada
nonlinearTseries::mutualInformation(dlserie, lag.max = 100,
n.partitions = 50,
units = "Bits",
do.plot = TRUE)
## $time.lag
## [1] 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
## [19] 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35
## [37] 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53
## [55] 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71
## [73] 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89
## [91] 90 91 92 93 94 95 96 97 98 99 100
##
## $mutual.information
## [1] 3.908374 1.421193 1.455035 1.402739 1.493103 1.555906 1.622355 1.429019
## [9] 1.455684 1.472721 1.553306 1.524770 1.539951 1.462847 1.493069 1.532704
## [17] 1.674945 1.637317 1.606369 1.709394 1.643386 1.664606 1.578555 1.629250
## [25] 1.658723 1.653978 1.706768 1.625568 1.645700 1.740113 1.544253 1.622424
## [33] 1.684184 1.678909 1.651387 1.633178 1.733924 1.629698 1.774529 1.713856
## [41] 1.833378 1.693334 1.762807 1.807361 1.655236 1.800483 1.853089 1.751038
## [49] 1.714328 1.654343 1.759697 1.794311 1.768562 1.839690 1.937591 1.927800
## [57] 1.721239 1.820051 1.912216 1.967309 1.800184 1.839622 1.729774 1.779251
## [65] 1.825002 1.785359 1.792801 1.772154 1.735934 2.024498 1.849860 1.820850
## [73] 1.911494 1.844747 1.830496 1.857667 1.804805 2.060829 1.956765 1.879505
## [81] 2.010459 2.004286 2.163420 1.962755 1.773337 1.823701 1.839891 2.023176
## [89] 2.049709 2.022580 2.166978 2.020244 2.327151 2.314320 2.222683 2.292481
## [97] 2.305321 2.231270 2.320423 2.370951 2.299530
##
## $units
## [1] "Bits"
##
## $n.partitions
## [1] 50
##
## attr(,"class")
## [1] "mutualInf"1.7 Exploración de la Estacionalidad
monthplot(dlserie)
tsibble_lserie %>%
na.omit()|>gg_subseries(diff_lserie,period=12) +
theme_minimal()
monthplot(dlserie)
tibble_lserie %>%na.omit()|>
mutate(
Mes = str_c("", as.character(lubridate::month(Fecha,label=TRUE)))
) %>%
plot_time_series(
.date_var = Fecha,
.value = diff_lserie,
.facet_vars = Mes,
.facet_ncol = 4,
.color_var = Mes,
.facet_scales = "fixed",
.interactive = FALSE,
.legend_show = FALSE,
.smooth = FALSE
)
ggseasonplot(dlserie)
1.7.1 Gráfico de cajas
tibble_lserie %>%
na.omit() %>%
plot_seasonal_diagnostics(.date_var = Fecha,.value = diff_lserie,
.feature_set = c("month.lbl"),.geom="boxplot")ggplot(tibble_lserie %>%
na.omit()|>
mutate(Mes = str_c("Mes ", as.character(lubridate::month(Fecha)))),
aes(x = diff_lserie)) +
geom_density(aes(fill = Mes)) +
ggtitle("Estimación de la densidad vía Kernel por mes") +
facet_grid(rows = vars(as.factor(Mes))) 
1.7.2 Periodograma
spectrum(as.numeric(dlserie))
Periodgramadlserie <- spectrum(as.numeric(dlserie),log='no')
ubicacionlserie=which.max(Periodgramadlserie$spec)
Periodgramadlserie <- spectrum(as.numeric(dlserie),log='no', main = "Periodogram")
ubicacionlserie <- which.max(Periodgramadlserie$spec)
abline(v = ubicacionlserie, , col = 'darkred', lty = 2)
sprintf("El valor de la frecuencia donde se máximiza el periodograma para la serie es: %s",Periodgramadlserie$freq[ubicacionlserie])## [1] "El valor de la frecuencia donde se máximiza el periodograma para la serie es: 0.391666666666667"sprintf("El periodo correspondiente es aproximadamente: %s",1/Periodgramadlserie$freq[ubicacionlserie])## [1] "El periodo correspondiente es aproximadamente: 2.5531914893617"1.7.3 Ajuste de la estacionalidad con componentes de Fourier y Dummy
tsibble_serie <- as_tsibble(serie)
diff_tsibble <- tsibble_serie|>
mutate(logdiff_serie = difference(log(value)))|>
select(logdiff_serie)
# Explore diferentes valores de K
Modelo_serie_diff<-diff_tsibble|>
model(Fourier1seriediff = ARIMA(logdiff_serie ~ fourier(K=2) +
pdq(0, 0, 0) + PDQ(0, 0, 0)))
real_ajustado1 <- diff_tsibble %>%
left_join(fitted(Modelo_serie_diff,by=index)) %>%
select(-.model)
real_ajustado1 %>%
autoplot() +
geom_line(data=real_ajustado1,
aes(y=logdiff_serie, colour="real"))+
geom_line(data=real_ajustado1,
aes(y=.fitted, colour="ajustado"))+
scale_color_manual(name = "real/ajustado",
values = c("real" = "black", "ajustado" = "red")) +
theme_minimal()
# Ajuste Dummy
Modelo_serie_diff_Dummy<-diff_tsibble|>model(
DummyAirdiff=ARIMA(logdiff_serie~season()+pdq(0, 0, 0) + PDQ(0, 0, 0))
)
Modelo_serie_diff_Dummy<-diff_tsibble%>%left_join(fitted(Modelo_serie_diff,by=index))%>%select(-.model)
Modelo_serie_diff_Dummy %>%
autoplot() +
geom_line(data=Modelo_serie_diff_Dummy,aes(y=logdiff_serie,colour="real"))+
geom_line(data=Modelo_serie_diff_Dummy,aes(y=.fitted,colour="ajustado"))+
scale_color_manual(name = "real/ajustado", values = c("real" = "black", "ajustado" = "red")) + theme_minimal()
# Varios modelos la mismo tiempo
ajuste_final_models<-diff_tsibble%>%model(
Fourier1Airdiff=ARIMA(logdiff_serie~fourier(K=1)+pdq(0, 0, 0) + PDQ(0, 0, 0)),
Fourier2Airdiff=ARIMA(logdiff_serie~fourier(K=2)+pdq(0, 0, 0) + PDQ(0, 0, 0)),
Fourier3Airdiff=ARIMA(logdiff_serie~fourier(K=3)+pdq(0, 0, 0) + PDQ(0, 0, 0)),
DummyAirdiff=ARIMA(logdiff_serie~season()+pdq(0, 0, 0) + PDQ(0, 0, 0))
)
glance(ajuste_final_models)FALSE # A tibble: 4 × 8
FALSE .model sigma2 log_lik AIC AICc BIC ar_roots ma_roots
FALSE <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <list> <list>
FALSE 1 Fourier1Airdiff 0.0429 19.0 -32.1 -31.9 -23.7 <cpl [0]> <cpl [0]>
FALSE 2 Fourier2Airdiff 0.0428 20.2 -30.4 -29.8 -16.4 <cpl [0]> <cpl [0]>
FALSE 3 Fourier3Airdiff 0.0413 23.4 -32.7 -31.7 -13.2 <cpl [0]> <cpl [0]>
FALSE 4 DummyAirdiff 0.0409 26.5 -29.0 -26.1 4.42 <cpl [0]> <cpl [0]>ajuste_final_models %>%
select(Fourier1Airdiff)%>%coef()FALSE # A tibble: 2 × 6
FALSE .model term estimate std.error statistic p.value
FALSE <chr> <chr> <dbl> <dbl> <dbl> <dbl>
FALSE 1 Fourier1Airdiff fourier(K = 1)C1_12 -0.0101 0.0268 -0.377 0.706
FALSE 2 Fourier1Airdiff fourier(K = 1)S1_12 0.0288 0.0266 1.08 0.282Modelo_serie_diff_models<-diff_tsibble%>%
left_join(fitted(ajuste_final_models)|>
group_by(.model)%>%
pivot_wider(names_from = .model, values_from = .fitted))
Modelo_serie_diff_models %>%
autoplot() +
geom_line(data=Modelo_serie_diff_models,aes(y=logdiff_serie,colour="real"))+
geom_line(data=Modelo_serie_diff_models,aes(y=Fourier1Airdiff,colour="ajustadoFourier1"))+
geom_line(data=Modelo_serie_diff_models,aes(y=Fourier2Airdiff,colour="ajustadoFourier2"))+
geom_line(data=Modelo_serie_diff_models,aes(y=Fourier3Airdiff,colour="ajustadoFourier3"))+
geom_line(data=Modelo_serie_diff_models,aes(y=DummyAirdiff,colour="ajustadoDummy")) +
scale_color_manual(name = "real/ajustado", values = c("real" = "black", "ajustadoFourier1" = "red","ajustadoFourier2" = "blue","ajustadoFourier3"="green","ajustadoDummy"="yellow"))